c=========================================================
C  Make the Green's function on X mesh
      subroutine gf1(ln,eproj,eta,ucentr,cntfug,gridx,nx,
     >   jdouble,njdouble,regcut,expcut,reg,beta,jstart,jstop,phase)

      dimension ucentr(nx),cntfug(nx),gridx(nx),jdouble(njdouble),
     >   reg(nx),beta(nx)
      complex phase
      complex*16 xx,eta1,zlmin,cfc(1),cgc(1),cfcp(1),
     >   cgcp(1),sig(1)

      call reglr(ln,eproj,eta,ucentr,cntfug,gridx,nx,
     >   jdouble,njdouble,regcut,expcut,reg,jstart,jstop,phase)

      xlmin=ln
      zlmin=cmplx(dfloat(xlmin),0d0)
      ifail=1
      nl=1
      
C  Now to find the iregular solution. We take care so as not to start
C  at values that are too small. Use routine COULCC to get the G function
      j=jstop
      if (eproj.lt.0.0) then
         phase=phase/(0.0,1.0)
         mode1=12
         kfn=0
         eta1=eta/(0.0,1.0)
         do while (j.ge.jstop-1)
            xx=sqrt(cmplx(dfloat(eproj)))*dfloat(gridx(j))
            call coulcc(xx,eta1,zlmin,nl,cfc,cgc,cfcp,cgcp,sig,mode1,
     >         kfn,ifail)
            beta(j)=abs(cgc(1))
            j=j-1
         end do 
         phase=phase*cgc(1)/abs(cgc(1))
         rph=real(phase)
         aph=aimag(phase)
         if (abs(aph).gt.1e-5) then
            print*,'Problems in making Green''s functions'
            print*,'Imaginary part of phase is too big, expect 0',aph
         end if 
      else if (eproj.eq.0.0) then
         do while (j.ge.jstop-1)
            if (eta.eq.0.0) then
               beta(j)=gridx(j)**(-ln)
            else
C  Here to evaluate the Coulomb G.F. at zero energy
               xx=dcmplx(sqrt(8.0*gridx(j)))
               mode1=2
               kfn=2
               xlmin=2*ln+1
               zlmin=cmplx(dfloat(xlmin),0d0)
               call coulcc(xx,eta1,zlmin,nl,cfc,cgc,cfcp,cgcp,sig,
     >            mode1,kfn,ifail)
               beta(j)=-real(cgc(1))*sqrt(gridx(j))
            end if
            j=j-1
         end do 
      else
         jj=jstop
         do while (j.ge.jstop-1.or.abs(beta(jj)).lt.0.1)
            xx=sqrt(cmplx(dfloat(eproj)))*dfloat(gridx(j))
            eta1=cmplx(dfloat(eta),0d0)
            mode1=2
            kfn=0
            call coulcc(xx,eta1,zlmin,nl,cfc,cgc,cfcp,cgcp,sig,mode1,
     >         kfn,ifail)
            beta(j)=real(cgc(1))
            jj=j
            j=j-1
         end do
      end if 
      
      j=j+2
      s2=beta(j-1)
      s3=beta(j)
      call nmrvb(ln,eproj,ucentr,cntfug,gridx,nx,
     >   jdouble,njdouble,s2,s3,beta,jstart,j)
      t2=1.0/real(phase)
      t=t2*imag(phase)
      if (abs(t).gt.1e-4) then
         DO J=JSTART,JSTOP
            BETA(J)=beta(j)*t2-t*reg(j)
         end do
      else
C  T2 may be -1.0         
         DO J=JSTART,JSTOP
            BETA(J)=beta(j)*t2
         end do
      end if 
      RETURN
      END
